Trends in Resonance Energy Shifts and Decay Rates for Bose Condensates 

in a Harmonic Trap 
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This is a study of quasi-discrete Bogoliubov quasi-particles in a spherically symmetric harmonic 
trap. We first evaluate analytically the aymptotic energy shifts of the high energy modes and 
find them to have l/s/n dependence on the number of radial nodes, n, consistent with earlier 
semiclassical discussions. To address the question of the widths or decay rates, 7, we attempt to 
clarify previous discussions by deriving an implicit equation for 7 from an assumption of exponential 
decay. Numerically, we study the trends in the behavior of the widths as a function of temperature, 
energy, particle number and scattering lengths. In particular, we find that the width due to Landau 
decay rises rapidly at low n and then declines, while the Beliaev decay rate rises slowly with n. As 
temperature — * 0, Beliaev decay reaches a constant (> for n > 0), while the Landau decay rate 
goes to zero. The decay rate, 7, is approximately linear in the s-wave scattering length. 



PACS numbers: 03.75.Fi, 05.30. Jp, 05.10.-a, 67.40. Db 



I. INTRODUCTION 

In many finite quantum systems, such as atoms, nuclei, 
and molecules, an extensive structure of excited states 
plays an important role in dynamical situations and for 
diagnostics. A homogeneous quantum fluid, such as liq- 
uid 4 He, exhibits a few excitation modes whose energy 
as a function of wavenumber and temperature has been 
intensively studied jlj. Harmonically trapped, dilute al- 
kali atom Bose condensates have recently been found 
to have at least two prominent excitation modes |^-||] 
and theoretically there is the full manifold of Bogoliubov 
quasi-particle modes. In calculating thermal averages, 
one can perform a sum over such discrete quasi-particle 
modes and obtain results comparable to the more com- 
monly used semiclassical integral over momentum space 
H . This raises the issue of the width of higher- lying res- 
onances, which is accentuated by recent observations of 
surface modes |7j and a chance observation of a longitudi- 
nal mode H. In other words, if experimental conditions 
are such that higher modes can be excited, under what 
circumstances are their widths narrow enough that they 
can be resolved? Could collective or single-particle ex- 
cited states of Bose condensates become research tools 
as useful as Rydberg states of atoms? 

We first focus on the energies of excited modes. Ana- 
lytic results have been obtained |j- 11 most notably from 
the hydrodynamic approximation, which involves various 
approximations to the kinetic energy. For highly excited 
states, these approximations break down. Without mak- 
ing the hydrodynamic approximation, Refs. jl2| jl5| used 
Bogoliubov theory to numerically compute low-lying ex- 
citations. High- lying levels were calculated numerically 
by You and Walsworth [ 16 ] . Semiclassical methods have 
enabled Csordas et al. [17, ITU to obtain precise expres- 
sions for the shifts (relative to bare harmonic oscillator 
states) of high lying excitations in a spherically symmet- 



ric trap. The more general case of cylindrical symmetry 
was found to lead to chaotic motion. In this paper, we 
present a very simple derivation of the lowest order ap- 
proximation to the energy shifts, somewhat in the spirit 
of the quantum defect parameterization of excited Ryd- 
berg levels of atoms, and we also present numerical re- 
sults from Bogoliubov-Popov theory. 

Since the first experimental observations of resonance 
excitation and decay in harmonic traps a num- 

ber of theoretical approaches have been presented on the 
question of resonance width or decay [^9|-^l|. Much of 
this effort has been directed toward the JILA data, al- 
though more recently Landau and also Beliaev damping 
of a "scissors" mode in a TOP trap has been observed 
at Oxford Arguably the most successful approach 

in explaining the JILA data has been kinetic theory with 
Gaussian Ansatz for the condensate and excited mode 
plus collision integrals for the condensate-noncondensate 
interaction p^psf . It was concluded that at a certain 
temperature the observed mode changed character from 
out-of-phase to in-phasc oscillation of the condensate 
and the thermal atoms, and with this interpretation the 
calculations matched the experimental data quite well. 
Good agreement has been offered also by calculations 
with a stochastic sampling of classical orbits correspond- 
ing to quasi-particle quantum states [pl[ and more re- 
cently with the dielectric formalism [p9|. In the present 
work, we will not consider again this experimental data 
but instead address the question of the widths of high- 
lying excitations in a spherically symmetric trap. Our ap- 
proach makes use of accurate (numerical) quasi-particle 
modes, and is valid when the perturbation of the system 
is small. This approach has the advantage of making 
it possible to consider explicitly the decay channels for 
Landau and Beliaev decay processes. We will present 
numerical results for decay widths as a function of atom 
number, interaction strength and temperature. 
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Expressions for Landau and Beliaev decay of reso- 
nance excitations of trapped Bose eases have been pre- 
sented by Pitaevskii and StringariJ 20 L by Guillemas and 
Pitaevskii |30|, and by Giorgini ]25| , [26| ]. Effectively, we 
adopt the expressions given in these papers, but pro- 
vide additional justification. Our discussion illuminates 
some of the many-body aspects and the analogy with 
Weisskopf-Wigner theory of radiative decay. One diffi- 
culty with previous expressions is that the terms in the 
sum over decay modes involved a delta function in energy, 
which, if interpreted exactly, will tend to be null for dis- 
crete quasi- particle states. By introducing the Weisskopf- 
Wigner Ansatz of an exponentially decaying initial state, 
we obtain decay terms with a Lorentz profile factor in 
place of the energy delta function. 

In this paper, our numerical results are intended not 
primarily to establish precise values for the energy shifts 
and the damping widths, but instead to seek patterns 
and trends as a function of various parameters. We con- 
sider only a spherically symmetric trapping potential. 
An outline of the paper is as follows. In Sec. II, we 
review the basic equations of Bogoliubov-Popov theory. 
In Sec. Ill, we use time independent perturbation the- 
ory to establish analytically, within the Thomas-Fermi 
approximation, the asymptotic behaviour of the energy 
shifts. In Sec. IV. A we review the basic theory of ex- 
ponential decay in time-dependent perturbation theory 
and then proceed to derive an implicit equation for the 
decay width. In Sec. IV. B, the lowest order interaction 
terms that lead to damping of an oscillating condensate 
are derived. Sec. IV. C uses these results in conjunction 
with the methods developed in Sec. IV. A to derive an 
expression for the decay rate of an oscillating condensate. 
Section V presents numerical methods and results of nu- 
merical computation for the width and its dependence on 
various physical parameters. 



II. HAMILTONIAN AND BOGOLIUBOV 
EQUATIONS 

The second-quantized grand-canonical Hamiltonian for 
a weakly interacting Bose gas in a spherical harmonic 
trap is 

K = J drft (H ho - n)$ + | / dr fiftM, (1) 
where fj, is the chemical potential and 

-h 2 v 2 1 



2m 



2 2 

-mtou^r 



g = Airh a/m. (2) 



with a the s-wave scattering length and Who/27r the har- 
monic frequency of the trap. The field operators are de- 
composed in the usual way into a macroscopic condensate 
wavefunction and a noncondensate operator 

4>(r, t) ~ $(r, t)) + $(r, t) = $(r, t) + 4>(r, t). (3) 



The temporal dependence of the condensate wavefunc- 
tion is attributed entirely to a small amplitude oscilla- 
tion (fluctuation) of the condensate about its equilibrium 
value <&o(r) 



$(r,t) = $ (r) + 6$(r,t). 



(4) 



These substitutions are made and the resulting expres- 
sions linearized with respect to the fluctuations. Sub- 
sequently equilibrium mean fields ny(r) = ($(/>) and 
mj-(r) = (4>4>) are introduced to reduce the Hamiltonian 
to a quadratic form in the noncondensate operators. If 
<f>o( r ) satisfies the Gross-Pitaevskii equation 



[H ho - /i + .g(|$ | 2 + 2n T )] $ = 0, 



(5) 



the first order terms in cf> cancel out. The present work 
does not deal with the regime near the critical tempera- 
ture, T c , and thus the Popov approximation will be ade- 
quate and so terms in mx are neglected. 

If the coupling between the fluctuations and the ther- 
mal operators are left out, the rest of the Hamiltonian is 
diagonalized by the Bogoliubov transformation: 

k*,t) = J2 [u k (t)a k (t) + u£(r)4(*)] 

k 

ft 0", *) = Yl ["* (*) + Vk (*)] ' ( 6 ) 

k 

provided that the functions u k (r) and v k (r) obey the Bo- 
goliubov equations, 

[H ho -/i + 2.g(|$ | 2 + n T )] u k (r) + ,g(|$ | 2 H(r) 

= e k u k {r) 

- [H ho -fi+ 2g(|$ | 2 + n T )) v k (r) - 5 (|$ | 2 H(r) 

= € k v k {r) (7) 

In spherical symmetry, quantum numbers n, £, m may be 
attached to the quasi-pa rticle sta tes k. In what follows, 
we will scale r by tq = y^W^ho and energies by Uui^q. 

Numerical methods using the discrete variable repre- 
sentation (DVR) |34j have been applied |35|,|36| to solve 
the Gross-Pitaevskii and Bogoliubov equations. A ther- 
mal sum over discrete quasi-particle states is supple- 
mented at high energies by an integral over a scmiclassi- 
cal representation. 

III. ASYMPTOTIC ENERGY SHIFTS 

We first consider the energy shifts of the collective and 
single particle modes of harmonically trapped conden- 
sate, relative to the energies of the harmonic oscillator 
modes. We are interested in the asymptotic behavior 
of the shifts as a function of the angular momentum, I, 
and the number of radial nodes n of the excited modes. 
In this section we consider only zero temperature, which 
implies that almost all the atoms are in the condensate. 
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FIG. 1. Hartree-Fock Bogoliubov Popov (HFBP) and har- 
monic oscillator (HO) energies for £=0, 10 and 20, relative 
to the trap minimum. The chemical potential is shown as a 
horizontal line. The scaled scattering length is a/ro = 0.0072 
and number of atoms N = 20, 000. 

By way of orientation, we first discuss overall trends 
of the energies e n i, or more precisely, E n i = e n i + A*- I n 
Fig. [I], we have plotted numerically computed eigenval- 
ues E n i, for T — and for I — 0, 10 and 20 on an energy 
scale relative to the bottom of the harmonic well. On this 
scale, the condensate lies at an energy equal to the chem- 
ical potential pt, shifted up from the harmonic zero point 
energy by the repulsive atom- atom interactions, and all 
excited states lie above the condensate. "Bare" harmonic 
oscillator energies E^ are shown for comparison. We will 
be interested in how the values E n i approach the E nl , but 
first it is relevant to point out certain trends at low n. 

For £ > 0, centrifugal distortion pushes the excited 
state wavefunction away from the condensate and thus 
reduces the shift for very small n. For larger n, the inner 
turning point moves inwards and the overlap with the 
condensate, and hence the energy shift, increases over a 
certain range of n. Then spreading of the excited state 
wavefunction to larger values of r results in a decreasing 
shift. These effects are fully discussed in Ref. [(L8) in the 
semiclassical context. 

A detailed and rigorous derivation of the energy shifts 
by WKB methods has been presented in |M from a con- 
sideration of classical particle motion in the presence of 
the condensate. Trajectories that do not penetrate the 
condensate have excitation energies that are simply the 
harmonic oscillator energies, 

Etc 



(2n + £ + 3/2). 



(8) 



For trajectories that do penetrate the condensate region, 
there is a shift: 



E„ 



(9) 



(10) 



that is found to be |lq| 

1 [4 M (2n + I + 3/2 - /i) - (I - 1/2) 2 ] 3 / 2 
n£_ 37r [2n + £ + 3/2- v} 2 

It is possible to obtain the n — > oo limit of the above 
expression by simple perturbation theory. For this pur- 
pose, we set nr(r) = as appropriate for T = 0. We 



define a 2 x 2 matrix Hamiltonian corresponding to the 
left hand sides of the coupled equations Eq. ([?]): 



Hi 



#ho-M 
-H ho +fi 



2 5 |$ | 2 5 |$ | 2 
-<7|$ | 2 -2 5 |$o| 2 



(11) 



The zeroth order eigenvector corresponding to eigenvalue 
E ni ~M is (VWi(r),0) T , where ip n i m (r) is an eigenstate 
of the spherical harmonic oscillator Hamiltonian, H^: 

V2n\ 



1pnlm(r, 



L n 2 (r )re * r - 



4tt 
(12) 



We have used standard notation for the spherical har- 

monics Y n i m {6,4>) and Laguerre polynomials L n 2 (r 2 ). 
Considering the second matrix operator in Eq. jTl|), the 
lowest order energy shift is: 



AE nl = 2g(^ n i m \ |$ | Ipnlm)- 



(13) 



An estimate of the behaviour of the shifts for large 
values of n may be obtained from the Thomas-Fermi ap- 
proximation for the ground-state density, 



l*ol 



2/z-r 2 
25 



(14) 



and the asymptotic form of the Laguerre polynomials 



L n 2 (r 2 ) ~ —^e 2T r 

\/7T 



^ +1 ^n 2 cos 



2V^r - (Z + 1)- 



Stirling's formula gives an approximate value 
n! 1 e' + 3 



1 



(15) 



(16) 



which, for large values of /, can be further simplified to 



(17) 



We thereby obtain the asymptotic first order shift in en- 
ergy 



AE nl 



f2f> 



dr{2n-r 2 )[l+ (-1)' cos(4V^r)] 



2 (2fi) 3 c+^)('+ 1 ) 
37t V n 



(18) 



to leading order in n. In the limit £ <C n — > oo, the shift 
becomes 



3 



SrU 



2 

3^ 



(19) 



which agrees exactly with the asymptotic value of Eq. 

To compare the asymptotic value of the shift obtained 
numerically with the value obtained from Eq. (|l0|), and 
also with values from Eqs. (|l^) through (jig), we plot 
these respective values multiplied by yfn over a range 
of values of n for I = 0, 5 and 10 in Fig. || (for 
TV = 50, 000, a sc /r = 0.0073). For £ = and also £ = 1, 
the simple perturbation theoretic results happen to be 
closer to the numerical results than Eq. (|l0|), but for all 
high values of £ and increasingly so, the WKB result is 
much closer than the perturbation theoretic result with 
asymptotic Laguerre functions. However, the WKB re- 
sult does tend to approach the asymptote more rapidly 
than the numerical results. It should also be noted that 
the use of the more accurate expression, Eq. (|l6| ) for Stir- 
ling's formula produces an improvement over the usual 
approximation. 
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FIG. 2. Energy shifts AE multiplied by yfn are plotted ver- 
sus n, the number of nodes in the radial wavefunction. The 
solid lines represent the energy shifts obtained from numerical 
solutions of the Bogoliubov equations ([?]). The long dashed 
lines above the solid lines represent the energy shifts calcu- 
lated analytically from the WKB theory given in Eq. (10). 
The medium dashed lines show perturbation theory results 
of Eq. (17), while for the short dashed lines the more exact 
Stirling's approximation in Eq. (15) is inserted in Eq. (17). 
In each case, the scaled scattering length is a/ro = 0.0072 and 
number of atoms N — 50, 000. 

We have also calculated first order energy shifts, Eq. 
( |l3| ) with exact harmonic oscillator wavefunctions rather 
than the asymptotic forms given above. The results are 
extremely close to the numerical results obtained with 



quasi-particle wavefunctions, from Eq. (Q). 

The fact that the asymptotic shift is independent of 
£ is demonstrated in Fig. |3| where we have plotted the 
energy shifts multiplied by yfn for I = Q, 5, 10 and 20. 
These shifts can be contrasted with the behaviour of Ry- 
dberg states of closed shell atoms, for which the asymp- 
totic shift is R6i/n 3 , where the dipole polarizability is 
6f ° l ~ 3ad/4:£ 5 , with a d the dipolar polarizability. Thus 
for atomic Rydberg states, the asymptotic shift exhibits 
a sharp £ dependence that is lacking here. One other con- 
trast with the theory of atomic Rydberg states is that the 
ion core can be considered effectively localized at the ori- 
gin, while for excitations of a Bose condensate, a spatial 
integral, as in Eq. (|l8|), is essential. 
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FIG. 3. Energy shifts A£=HFBP(n, I) - UO(n,£) multi- 
plied by the yfn are plotted for different values of £, showing 
that the asymptotic shift is independent of £, and approxi- 
mately equal to 2yJ (2/i) 3 /37r. The scaled scattering length is 
a/r = 0.0072 and number of atoms TV = 50, 000. 

In cylindrical symmetry, WKB wavefunctions exhibit 



chaotic behaviour 
out as for sphcrica. 



37fl , and the analysis cannot be carried 
symmetry. The asymptotic shift may 
be obtained from perturbation theory in a manner analo- 
gous to the above discussion, and results in the following 
expression: 



9 - 

\x e 



■■y/u p n p ijj z n z 



(20) 



where lo„ and oj z are trap frequencies, and m is an exact 
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quantum number. However, n p and n z are exact only 
when there is no condensate present. 



IV. LANDAU AND BELIAEV DAMPING OF 
EXCITATIONS IN A TRAPPED BOSE GAS 

We now consider the damping of these collective and 
single particle excitations. An estimate of the width of 
the modes is essential to establish the feasibility of re- 
solving high-lying excitations of a Bose condensate. The 
primary source of damping is the interaction of the collec- 
tive excitation of the condensate with the noncondensate 
or thermal part. For high-lying excitations, the possi- 
bilities for both Landau and Beliaev decay increase with 
temperature, so it is interesting to obtain these rates ex- 
plicitly. 

A detailed exposition of the theory of Landau and Be- 
liaev damping of quasi-particle excitations for Bose con- 
densates has been presented by Giorgini Q . The expres- 
sions for Landau damping were in agreement with expres- 
sions given earlier by Pitaevskii and Stringari |20(| , but 
the derivation was considerably more detailed. Giorgini's 
derivation is applicable to Bose atoms in harmonic trap 
potentials as well as to homogeneous Bose gases. Our dis- 
cussion takes another point of view that results in expres- 
sions with a Lorentz profile. Consequently, all interacting 
levels within the linewidth of a resonant mode contribute 
to damping by Landau or Beliaev mechanisms. We note 
again that references and (2^] also present theories 
of Landau damping for trapped Bose atoms. 



mt)) = u(t,o)\m)- 



(24) 



From Eq. ([23|) it follows that the operator satisfies the 
differential equation 



ih^Lu(t,o) = v(t)u(t,o) 



(25) 



which can be integrated with the initial condition 
U(0, 0) = 1 to give 



U(t,0) = 1 



ft .10 



dt'V(t')U(t',0). 



(26) 



The final state can be expanded in a basis of stationary 
eigenstates of Hq 



|*(t)>=£a»(t)|*«>, 



(27) 



with coefficients 

o«(t) - (* n |tf(t,Q)|*(0)) 

i f* 

= <ig#(0))-^~ dt'{^ n \V(t')\m m )a m (t'). (28) 

m ^° 

At this point we make the simplifying assumption that 
the initial state is an eigenstate of Ho, 1^(0)) = 
Moreover, if the time scales are short enough and there 
are a large number of possible final eigenstates, it can 
be assumed to a very good approximation that a,i(t) ^ 
a n ^i{t), so that 



A. Damping from Perturbation Theory 

We first review some of the methods of time-dependent 
perturbation theory |3§| , |39| ] that will be used in our ap- 
proach to the calculation of the widths. It will be con- 
venient to work in the interaction picture, in which the 
time dependence of the operators and the states is given 
by 



6{t) = e iHat/hQ se -iHot/H 



(21) 



(22) 



where Os signifies a time-independent operator in the 
Schrodinger picture. The Hamiltonian is H = Hq + V. 
The differential equations that determine the time devel- 
opment are 



4™> 



-v(t)\nt)) 

[O(t),H ]. 



(23) 



The time development of the state from time t = to a 
later time t can be effected by means of a time-evolution 
operator U (t, 0) acting on the initial state 



(*) 



h jo 



dt'^Vit'^CLiit') 



(29) 



If the time scales are sufficiently short that changes in 
the initial state can be ignored, we can set cii(t) « 1, and 
Eq. (E9) can be solved for each n ^ i to obtain: 



\a n (t)\ 2 = 2\{H> n \Vs\y<) 



1 - COs{(tJ n - LOi)t} 
[h(uj n - LUi)} 2 



(30) 



We have reverted to the Schrodinger picture in order to 
make the time-dependence explicit. E n — huj n is the 
energy associated with the nth eigenstate of the unper- 
turbed Hamiltonian Hq. Then the usual procedure is to 
assume that t ^ h/(E n — Ei), so that even though t is 
still small compared with the inverse decay rate, one is 
justified in writing the quantity in brackets as a delta 
function. In this way, one obtains Fermi's Golden Rule: 



r,. 



dt 



2n 



\^> n \Vs\^ l )\ 2 8{E n -E t ). (31) 



The sum of the partial widths, J2 n ^ n = ^ 7 gives the 
total decay rate at t = 0. 

This expression serves for the evaluation of transition 
probabilities in a wide range of physical problems. How- 
ever, for application to decay rates, it requires that there 
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be an intermediate range of times that are large com- 
pared with the oscillation periods but small compared 
with r _1 . For a typical radiative decay process, for ex- 
ample a D line of sodium with T/uo w 2 x 10~ 8 , this is 
no problem. For the resonance excitations of Bose con- 
densates reported in Ref. however, IV 'u ~ 0.01 to 
0.1, which makes the 6 function in Eq. (p^) problemati- 
cal. Hence it is appropriate to introduce the method of 
Weisskopf and Wigner ||j|[l(|, which removes the upper 
limit on the range of t by dropping the assumption that 
Oj(t) ~ 1. This is done by introducing an exponential 
decay of the initial state: 



a 



(*) 



(32) 



Eq. ( |29| ) then becomes 
r-l 



a n (t) ~ -1 / dt'i^Vsme^-^'e 



1 _ e i(w„-w ( )t e -|rt 



trt' 



o 



h{(u n -L0i)+iT/2} 



(33) 



Conservation of probability for a complete set of states 
for a closed system 



i-k(*)i 2 = £i<a*)f 



(34) 



then gives 



„-rt 



E 



^{K l -^) 2 + (r/2n 

x [l + e" rt -2e" rt/2 cos{(cj m -w l )0 . (35) 

As it stands, this expression does not have any assump- 
tions about the time-scales involved. We will take the 
long term limit, t — > oo, in order to remove the explicit 
dependence on time 



r/2 = E 



\(*m\V S \9i 



r/2 



K - w,;) 2 + (r/2y 



(36) 



We have multiplied both sides by r/2 to obtain an 
implicit equation for the width that can be solved itera- 
tively with a knowledge of the stationary states and their 
energies. This is the general form we will use to obtain 
decay rates for excitations in a Bose gas. 



B. The Damping Interaction 

We now consider the operators and the basis states 
that are relevant for a Bose-condensed system. Since we 
will work within the grand canonical ensemble, Hq in the 
previous section will be replaced by Kq , which is the part 
of the grand canonical Hamiltonian in Eq. (^) that is di- 
agonalized by the Bogoliubov transformation, Eq. <M). 



The eigenfunctions are the quasi-particle basis states as 
given in equations Eqs. (||)and (]?]). The quasi-particles 
themselves are considered to be in a thermal distribu- 
tion determined by Kq, so that the relevant statistical 
operator is 



-PK 



Tr e-P^ 



(37) 



Eigenstates of Kq may be written in a Fock represen- 
tation of quasi-particle occupation numbers |tii,71 2 , • • •), 
and satisfy 



K \ni,n 2 , 



k 



CfcTlfc 



|ni,n 2 , 



(38) 



where e k = huj k is the energy and n k is the occupation 
number of the fcth quasi-particle mode. The creation 
and destruction operators for the quasi-particle modes 
arc defined by 



a k \ni, ■ ■ ■ ,n k , ■ ■ •) = \fnk~\rii, ■ ■ ■ , n k - 1, ■ ■ •) 

«lK, •■■)«*)■••) = vWTT|ni, • • • ,n k + 1, • • •). (39) 

They obey the standard Bose commutator relations and 
their time-dependence in the interaction picture follows 
fromEq. (pi 



a k (t) = a k e 



-iujkt 



4w 



a k e 



(40) 



In this representation, the matrix elements of the statis- 
tical operator are 



(rai,n2, • • •! e pK ° |ni,n 2 , 



2{ mi } (mi,m 2 ,---|e & K ° |mi, m 2 , • • •) 



IL 



n*E, 



-f3t k m k 



(41) 



We next define the interaction that causes damping. 
We start with the interaction term in Eq. (1), sepa- 
rate the condensate wavefunction as done in Eq. (|), 
introduce the small amplitude fluctuation of Eq. (|) and 
linearize with respect to the fluctuation to get 
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drip' if)ip w g 



'/*{ 



+$o$o[(2<5$ + 5&)ft + (5$ + 25&)<f>] 

+(6$ + 5&)ft(fHj)\ (42) 

Here, and in subsequent expressions, for the sake of 
brevity, the spatial and temporal dependencies are not 
shown explicitly where they are obvious. In this paper, 
we assume that the condensate fluctuation mode is iden- 
tical with one of the quasi-particle modes Eq. (^) , which 
we designate by the subscript 'os ' : 



G 



5$(r,t) = u os (r)a os {t) + v* os {r)&t s (t) . 



(43) 



We ignore the dynamics of the noncondensate, i.e. the 
noncondensate quasi-particles are assumed to be in static 
thermal equilibrium, with a population given by the Bose 
distribution function. As in fl30[ , we regard the effect of 
the oscillation as an increase of the quasi-particle popu- 
lation of the mode labeled 'os ' from its value, f® s in a 
thermal distribution to a macroscopic value, which we 
call n os . Furthermore, we assume that the amplitude of 
the oscillation is small enough that the thermal equilib- 
rium of the noncondensate particles is not affected. 

Before making the quasi-particle substitutions, we note 
that the terms that are zeroth order in the noncondensate 
operators cannot contribute to the decay of the fluctua- 
tion due to conservation of energy. The first order terms 
will also be constrained by conservation of energy so that 
each quasi-particle removed from the condensate fluctu- 
ation will be compensated by an identical quasi-particle 
added to the noncondensate; this essentially leaves the 
number of quasi-particles in the condensate fluctuation 
unchanged. Thus the lowest order contribution to the de- 
cay of the condensate fluctuation comes from the terms 
bilinear in the noncondensate operators. Hence, the 
damping interaction is defined by 



V{t) = gJdr<P a [5&^+2(6$- 



(44) 



We transform to quasi-particles and retain all the 
terms that do not violate conservation of energy. These 
are grouped into terms (labeled V + (t)) that create a 
quasi-particle in the condensate fluctuation and those 
that annihilate one (labeled V~(t)): 



V(t)=V+(t) + V-(t). 



(45) 



These terms include two processes: (i) Landau decay 
in which a quasi-particle from the condensate fluctuation 
interacts with a noncondensate quasi-particle to create 
a higher energy quasi-particle and (ii) Beliaev decay in 
which a condensate fluctuation quasi-particle decays into 
two lower energy quasi-particles. The two processes are 
schematically represented below 



E 



Ej Ei 



Landau 
Beliaev 



(46) 



V (t) correspond to these processes going from left to 
right while those included under V + (t) proceed in the 
reverse direction. 

At this point, we make the time-dependence of the 
quasi-particle operators explicit. Also, since we are con- 
sidering a stationary condensate the quasi-particle wave 
functions Uk and Vk are taken to be real. After some al- 
gebra and rearrangements of the indices under the sum- 
mation, we obtain 



id 

i,3 



3 « 



(47) 



with ujj = uj os — (uji ± ujj) and 

Lij — 2g J dr $ [uos(uiVj + u { Uj + V t Vj)+ 

VosiUtVj +u i u j +v i v j )\ 
Bij = g J dr $ [^(w^ + v^j + VjUj 

+v os (v j u i + + v^)] . (48) 



C. Calculation of the Width 

The expressions for the interaction derived in the pre- 
vious section can now be applied to the methods devel- 
oped in Sec. IV. A to obtain the decay rate. Since we 
consider finite temperatures, the system has to be de- 
scribed by a density matrix. We incorporate our consid- 
erations from the previous section into the definition of 
the density matrix. The condensate oscillation is treated 
as an increase in population of a particular mode, labeled 
with subscript 'os', to a macroscopic value n os 3> 1. At 
the same time the enhanced population n os is considered 
sufficiently small to leave the remaining modes well ap- 
proximated by a thermal distribution. Thus we define 
the initial density matrix to be 



P(0) 



E 



W {nk} \n 1 ,n 2 ..n os ..)(n 1 ,n 2 ..n os ..\. (49) 



It is characterized by a well-defined and fixed value of n OSl 
for the condensate oscillation mode, whereas the other 
modes are thermally distributed over all possible values. 
This is reflected in the statistical weight 



n 



k^os ' 



-/3e k n k 



(50) 



The time evolution of the diagonal matrix elements of 
the density matrix is given by 

(mi, ?7i2, ..|p(t)|mi,m 2 , ..) 

= mi, 7712, •• ,"7-2, ••) 

= W i - n ^\(m ll m 2 ,..\U(t,0)\n 1 ,n 2 ..n os ..}\ 2 , (51) 

k^os 

where, U(t, 0) is the time evolution operator defined ear- 
lier. We will use it in the linear approximation 



tm,o)~i--i / dt'v(t'). 

h Jo 



(52) 
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The states |mi, m 2 , • • •) that evolve from the initial state 
can have arbitrary occupation numbers for all modes. For 
damping of the condensate oscillation mode, however, 
only states with m os ^ n os are relevant, for which case 



The sum over {n^} with the statistical weights essen- 
tially reproduces a thermal distribution of the modes of 
other than the condensate oscillation mode. In order to 
evaluate this, we consider a generic term, 



(m 1 ,m 2 , ..\p(t)\mi,m 2 , 



k^os 



dt' (mi , m 2 , .. | V(t') \m, n 2 ..n os ..) 



(53) 



£ 

ni,n 2 ,' 



n* 



-/3e k n k 



ill 



in, 



-/3c, mi 3 



Now we introduce an Ansatz similar to the Weisskopf- 
Wigner hypothesis that the number of particles in the 
condensate oscillation mode decays exponentially, i.e. 

n os (t) = n os e~ rt (54) 

This it Ansatz and the interaction derived in Eqs.( [45|) - 
( p8| ) allow us to evaluate the matrix elements in Eq. ( |53| ) 

ft „ 2 

' dt' (m 1 ,m 2 , ■ ■ ■ \V + (t')\n ll n 2 ..n os ..) = 



IL Em, 1 



-/3e,m, 



y 



rije 



l 



(58) 



Thus the thermal sum has the effect of replacing 

n -> fj = [e^ - (59) 

We divide by n os to obtain the probability of the 
change per initial quasi-particle: 



Li 



dt'e ( 



rit' 



t 2 



dt'e { 



(n os + l)(nj + l)ni 
(n os + XjUjUi 



7i os — n os + l 

m 3 — n j + 1 

rrii—ni — 1 

-n os —n os -\-l 
nij—nj — l 
rrii—ni — 1 



P(M;-1) 



*2 £ 



Lf 



(/i - fi) 



(60) 



(55) 



+5 



^ /e( ia,+ -ir)t 



(1 + fi + fi) 



t ^ 2 

dt'(mi,m 2 ,- ■ ■ \V~~ (t')\n 1 ,n 2 ..n os ..) 



Li 



-Bt 



^(k-y-jiy 



n os nj(rii + 1) 



The time integrals can be treated in the same way as 
in our general discussion in Section IV A . In the long 
time limit i.e. t — > oo, these integrals tend to 



nij—rtj — l 
rai—ni+1 



dt'e 



i(iu±.-lT)t' 



1 



(a,±) 2 + (r/2)2- 



(61) 



±r)t' 



«os(%- + + 1) 



m os — n os —1 
mi— ni + l 



(56) 



It is obvious that for any given set of occupation numbers 
{mi, m2..}, only one of the above terms will differ from 
zero. 

We add up all the non-zero elements that arise from 
Eq. ( p6| ) considering all possible values of m\,m 2 , 
and subtract from them all the non-zero elements that 
arise from Eq. ( pq ) . This gives us the probability for the 
number of particles in the condensate oscillation mode 
decreasing by one after time 'f when we start with 77. os 
particles in that mode. We use the macroscopic nature of 
the oscillation mode, which means n os is large compared 
to 1, to set 7i os ~ n os + 1 and obtain 



As for Eq. (ft6|), in the limit t — ► oo, the condensate oscil- 
lation has damped out completely so that P(l, t; —1) — > 
1, and we can obtain an implicit equation for T. At this 
point, we introduce 7 = T/2, which corresponds to the 
width parameter 7 used in references pQE5pO|. Thus 



1 



*£ 

hi 



Li 



7 



K) 2 +7 : 



(fi - fi) 



-Bi 



7 



(1 + fj + h) 



(62) 



P(rW;-l) 



k^os 



+Bl 



Li 



dt'e { ' 



1,1 (iuj. - — A 

dt e y 2 



r)t 



(n.j - Hi) 



where we have multiplied both sides by 7 = T/2 as in Eq. 
(|36|). In the numerical work, we will use this expression 
iteratively to calculate the width. 

To deduce an approximate relation for the dependence 
of 7 on the coupling parameter, g, and hence on a, the 
scattering length, we consider just those values for 
which tofj < 7. For this subset of interactions, which 
include the largest contributions to the width, we have 



(1 + rij + m) 



(57) 



hi 



[L% {f s - + B% (1 + fj + f^] . (63) 
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Since Ljj and Bij from Eq. (|4q ) are each linearly depen- 
dent on g, we conclude that 7 is approximately linearly 
dependent on g and a. 



D. Comparisons and Caveats 

Previous expressions for the damping rate given in 
Refs. [ pO| , p5|j3C| employed a delta function, as in Fermi's 
Golden Rule, Eq. ( 3l|). For example, for Landau damp- 
ing, Eq. (41) of Ref. M is 



1L 



(64) 



For a manifold of discrete states, this delta function 
would never be precisely satisfied. By assuming exponen- 
tial decay, we have obtained a width expression involving 
a Lorentz function with a finite width, that can be solved 
iteratively for 7. Our discussion in Sec. IV. A also makes 
clear that the delta function is difficult to justify when 
the ratio 7/w approaches values of 0.01 to 0.1 as is the 
case for excitations of Bose condensates. However, one 
can formally obtain a similar expression from Eq. (J64|) by 
replacing ujq, the unperturbed resonance frequency, with 
the perturbed lj, which is then assumed to be distributed 
over a range of values characterized by a Lorentz profile 
with a width 7. 

In the above derivation, we have explicitly calculated 
the decay of diagonal density matrix elements. However, 
if we consider d os as a c-number, a os = aj s = V^os, and 
evaluate the time evolution of the squared amplitude of 
the condensate wavefunction 

|$(r,t)| 2 H$ (r) + ( 5$(r,<)| 2 = 
I $ (r) 1 2 +2$ (r) ^NZ[u os (r)+v os (r)] cos(uj os t)e-^ 
+N os [u os (r) 2 + v os (r) 2 

+2u os (r)v os (r) cos(2ui os t))e~ rt , (65) 

we see that the "diagonal" term that decays as e _rt is 
second-order in the small amplitude fluctuations. The 
term linear in 3>o( r ) and linear in the fluctuations must 
be considered the primary observable oscillating function. 
This term oscillates with a frequency uj os and decays as 
e~ 7 * . It can be considered to arise from a coherent super- 
position of the condensate and the condensate fluctuation 
mode. For our derivation to be relevant to this term, one 
must assume that indeed 7 of Eq. ( |65| ) is equal to 7 of 
Eq. d62|), i.e. that the "off-diagonal" element decays at 
one-half the rate of the "diagonal" element. 

One can question whether the decay always follows a 
simple exponential law as we have assumed. For a more 
precise estimate of the time evolution of the resonance 
oscillation, numerical integration of coupled equations, 
including important off-diagonal density matrix elements 
such as the coherent term in Eq. (pq), would be needed. 



We note also that we have neglected the thermal atoms 
in the os mode, which will be significant at certain tem- 
peratures. An accurate calculation of the contribution 
of this mode to the damping process raises questions we 
have not resolved. However, in view of the additional 
constraint imposed by the Lorentz profile, we estimate 
that Landau damping by the os mode itself is not signif- 
icant. 



V. NUMERICAL RESULTS FOR EXCITATION 
WIDTHS 

To obtain the energies e = Tilj and the quasi-particle 
wave- functions Uk and the Vk, we solve self-consistently 
the Gross-Pitaevskii equation for the condensate and the 
Bogoliubov equations, as given in Eqs. (Eh and (Q). As 
stated above, here we accept the Popov approximation 
and neglect the anomalous density of the thermal par- 
ticles, but always nr(r) is included. These equations 
are solved numerically using the discrete variable repre- 
sentation (DVR) p4p5| ], a method derived from Gaus- 
sian quadrature methods. The sums over discrete quasi- 
particle states are supplemented by local density approx- 
imation (LDA) expressions which treat the high- 
lying states by a semiclassical integral over momentum. 
More details of this part of the calculations can be found 
in (36 1 . 

To compute the widths we put in a trial value of 7 on 
the right side of the implicit Eq. (^) and evaluate the 
new value of 7. The calculation is iterated until the input 
7 is equal to the calculated 7 within a pre-defined toler- 
ance. We made careful tests to ensure that the resultant 
widths are independent of various computational param- 
eters. In particular, in evaluating the sum over states 
we include all states such that Lof^ fell within a specified 
range. For each mode studied, the range was increased 
until the width was insensitive (within the tolerance) to 
the increments. Also we verified that the width was in- 
sensitive to initial trial values. Most of our calculations 
were done for i = excitations. The set of perturbing 
levels ranged from I — up to i max , ranging from 20 
to 80, which was found sufficient to obtain a converged 
value of 7. 

The above computational procedure is in principle 
more precise than the methods used in Ref. |30|, which 
employ condensate and excited state wavefunctions com- 
puted at T = with populations appropriate to the 
stated T. Although therefore we do not expect exact 
agreement with these results, we made extensive com- 
parisons to provide qualitative tests of our approach. For 
this purpose we used the same atom number values and 



scaled scattering length, a/r = 7.36 x 10 3 as in Q 
For example, the temperature variation of the widths for 
N = 50,000 and 150,000 from [|o) and from our results 
is shown in Fig. 0. 
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FIG. 4. The decay width, 7, for the lowest energy mode 
(n = 1), computed in this work as function of tempera- 
ture is compared to the widths computed by Guilleumas and 
Pitaevskii (GP) in Ref. [30]. The width is scaled by the 
monopole resonance frequency S7ho — v5who an d the tem- 
perature is scaled by the chemical potential, /1, as in this 
reference. 

From our numerical calculations for a spherically sym- 
metric harmonic trap we observe the following trends: 

(i) As in the data presented in ||(| , we find that the to- 
tal width typically increases with temperature somewhat 
linearly, as shown in Fig. |^. For the higher resonances, 
non-zero values at T=0 reflect the Beliaev damping rate. 




kT/( ttc V) 

FIG. 5. The decay width 7 as function of scaled tempera- 
ture, is shown for several condensate modes, for N=150,000 
and a/r = 0.00736. 

(ii) Figure ^ explains that the temperature increase of 
the width is due especially to the Landau contribution, 
while the Beliaev part increases much more slowly. This 
is plausible since the Landau process involves the inter- 
action of a thermally populated quasi-particle mode with 



the oscillating mode in question, while the Beliaev pro- 
cess is simply a decay into two modes of lower energy. At 
higher temperature, of course there are more thermally 
populated modes available. 




kT/( ttc °ho) 

FIG. 6. The total decay width and the Landau and Beliaev 
contributions as function of scaled temperature kT/hujho are 
shown for the quasi-particle mode with n = 15, iV=150,000 
and a/r = 0.00736. 
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FIG. 7. The decay width 7 as function of the number of ra- 
dial nodes, n, shown for a wide range of scaled temperatures, 
for iV=150,000 and a/r = 0.00736. 

(iii) The total width rises sharply for low n, the radial 
quantum number, and then reaches a nearly flat plateau 
at a value of n roughly proportional to T, as shown in Fig. 
[?]. We have correlated this effect with the factor /? — /? 
in Eq. (fS2|). For the (j, i) values that contribute most 
to 7, at higher temperatures, /? is more nearly equal to 
ff. The temperature scale at which this occurs is set 
approximately by {Ei — Ej)/k. At higher n values, is 
negligible compared with /?, so the width as a function 
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of n flattens out. Note in this figure that for low values of 
the temperature, the width is small compared with fwJho 
even for high values of n. 



However, if means can be found to achieve a perfectly 
spherically symmetric harmonic trap and then to excite 
a singlet I manifold, our conclusion is that the highly 
excited levels would be resolvable in this ideal case. 




20 30 



FIG. 8. The total decay width and the Landau and Beliaev 
contributions as function of the number of radial nodes n of 
the condensate fluctuation are shown for scaled temperature 
kT/hu; ho = 12, iV=150,000 and a/r = 0.00736. 



(viii) Although the above results were entirely for the 
case I = 0, we have performed calculations for I up to 6. 
For the cases studied, have found only a slight variation 
of the width with i. 



0.16 
0.12 



a/r n = 0.00736. 
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(iv) Figure |g, again breaks this trend down into Landau 
and Beliaev components for one particular case. The ap- 
parent plateau occurs because of a near balance between 
a decrease of the Landau rate and an increase in the Be- 
liaev rate. Since for higher energy, there are more final 
states available, one does expect an increase in the Beli- 
aev rate, but it is surprising that the increase is as grad- 
ual as it is. Also, the n dependence of the Landau rate 
is rather unexpected. We attribute the decrease at large 
n to the decrease of the overlap integrals in Eq. (48), 
since the excited state wavefunction will be increasingly 
extended in space. 

(v) The total width increases with scattering length 
in approximately linear fashion, as shown for one set of 
conditions in Fig. @. This is consistent with Eq. (63). 
A larger scattering length implies stronger interaction 
among the particles and hence more rapid damping. 

(vi) In Fig. |l0|, we have plotted widths for three differ- 
ent particle numbers. With fewer particles, the conden- 
sate is smaller and also the quasi-particle wavefunctions 
are closer to bare harmonic oscillator functions. For both 
reasons, the integrals in Eqs. (^) yield smaller values for 
the coupling parameters. 

(vii) Especially for relatively low temperature and 
smaller values of the scaled scattering length, the widths 
are relatively small compared to the energy level separa- 
tion, even at high values of energy or n. Unfortunately, if 
one attempts to excite high energy levels, typically many 
different I values would be excited and would overlap. 



FIG. 9. The decay width 7 as function of number of radial 
nodes, n, of the condensate fluctuation, shown for two differ- 
ent values of the scaled scattering length a/ro. The data are 
for N = 50,000, and kT/hu^c — 16. Also shown is the ratio 
of the widths (right axis), which may be compared with the 
value of 6.13 for the ratio of the scattering lengths. 
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FIG. 10. The decay width, 7, as function of number of 
radial nodes, n, of the condensate mode, is shown for three 
different values of the number of particles, N. The data are for 
scaled scattering length a/ro — 0.00736 and kT/huJho = 12. 



11 



VI. DISCUSSION AND CONCLUSIONS 

A key question motivating this study has been the ex- 
tent of validity of the quasi-particle representation for 
trapped Bose atoms. In light of recent kinetic theory 
results H^j2^], it may be that the quasi-particle repre- 
sentation is not essential for considering resonance de- 
cay. The appealing but complex picture of Landau and 
Beliaev decay processes may not be necessary. However, 
the Gross-Pitaevskii equation and the Bogoliubov quasi- 
particle approach to Bose condensates has been shown 
to be a powerful method for modeling thermal equilib- 
rium situations. Also in view of past successes pl| , |l5| ], 
it seems worthy of further study as well as application 
to more highly excited resonances as we have done here. 
The dielectric formalism ]2j| employs quasi-particle basis 
states, but goes beyond the present approximations. We 
are currently analyzing these other approaches to make 
comparisons. 

It is useful to distinguish between discrete quasi- 
particle modes, characterized by eigenvalues and 
eigenf unctions Ufc(r) and Wfc(r) from Eq. (Q), and semi- 
classical quasi-particle functions, u(p, r), v(p, r) as used 
in |J. The widths that we have computed for reso- 
nantly excited states may be presumed to apply also to 
the thermal quasi-particle states, which must thus be 
considered quasi-discrete quasi-particle (QDQP) states. 
When the widths of the QDQP states are comparable 
to their spacing, one can conclude that discrete basis 
states are not relevant, and a local semiclassical represen- 
tation will suffice. Spherical symmetry presents a special 
case because the single-particle excitations at sufficiently 
high energy within a given I manifold have a spacing of 
2hujho- We do find that for TV = 150, 000 for example, for 
T < 12hidho/k, typical widths are < O.lOfiWho, even for 
high- lying states with n up to at least 30. Going from 
spherical to cylindrical symmetry, the spacing between 
levels within a given m manifold decreases with energy, 
and thus semiclassical quasi-particle functions are more 
likely to be accurate. Consistent with this, we have ob- 
served that the semiclassical quasi-particle functions be- 
come more reliable at lower energies in cylindrical than 
in spherical symmetry for calculating thermal sums over 
quasi-particle states p6| , p2[ . 

Very recently, a harmonic trap close to spherical sym- 
metry has been produced |43| ], and we are informed ]44| ] 
that axial and transverse frequencies can be made equal 
within about 1%. Although we have found above that 
resonances within a given £ manifold are under certain 
circumstances resolvable, to our knowledge no one has 
yet devised an excitation method that would be I selec- 
tive, as are electric dipole excitations of atomic Rydberg 
levels. One can hope that such techniques for this pur- 
pose will be developed. You and Walsworth |l(| have 
suggested that time-varying magnetic fields be config- 
ured so as to excite specific modes. With the help of 
focused laser beams, the MIT group has observed reso- 



nances with n z as high as 8 in a cigar-shaped trap, and 
surface modes |8| were excited by a sharply focused per- 
turbing laser beam, moved rapidly around the transverse 
surface of the "cigar" trap. The technique of stimulated 
Bragg scattering |45|-|47| offers the capability to select 
the excitation energy and wavenumber, but is not £ se- 
lective in its current form. Possibly by special tailoring 
of the laser beams, one could achieve a nearly spherically 
symmetric excitation function. 

Our review of the underlying theory of Bose conden- 
sate damping processes has exposed some interesting 
contrasts with the Weisskopf-Wigner theory of radiative 
damping. The decay rates of Bose condensate excita- 
tions are typically not a small fraction of the frequency, 
which makes the usual delta function limit problemati- 
cal. Furthermore, the decay of interest is very likely the 
"off-diagonal" density matrix element between the con- 
densate and the excited mode rather than the excitation 
intensity. In the language of NMR, we have calculated 
here Ti, but T% may be what is observed. To model 
these various circumstances in a more accurate manner, 
it may be necessary to explicitly integrate coupled time- 
evolution equations for various diagonal and off-diagonal 
elements. Also, the assumption of simple exponential 
decay may not always be valid. Modulated decays have 
been reported in calculations with an assumed ID geome- 
try [Q and also from Monte Carlo procedures [fl9) . As is 
well known, a system composed only of discrete frequen- 
cies will exhibit recurrences rather than damping. This 
conclusion may not apply if the number of interacting 
levels is effectively infinite, but at sufficiently low tem- 
perature, there may not be effectively an infinite number 
of accessible quasi-particle levels. Again, more detailed 
modeling of the time evolution may be needed. 

In summary, we have presented numerical and pertur- 
bation theory results for excitation frequencies of Bose 
particles in a spherically symmetric harmonic trap, and 
compared with previously derived expressions. The ener- 
gies approach the bare harmonic oscillator values asymp- 
totically as (4/37r)/i 3 / 2 / \/2/n independent of£ (for n >> 
£), where n is the number of nodes in the radial equa- 
tion. The asymptotic regime is reached more slowly for 
higher £. Secondly, we have presented a derivation of 
the excitation widths in agreement with previous discus- 
sions, but based from the beginning on an assumption 
of exponential decay, with an ensemble of discrete quasi- 
particle states, hence yielding a Lorentzian width expres- 
sion rather than a delta function. Thirdly, we have pre- 
sented numerical results for the widths over a range of 
parameter values. They showing that the Landau pro- 
cess rises more rapidly with temperature, while the Be- 
liaev process rises more rapidly with n, except for very 
low n. Also, compared with the frequency separation of 
2Hlo within an I manifold, the widths are often small, 
especially for low temperatures and smaller values of the 
scaled scattering length. 
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